*This program runs a probit and calculates the associated marginal effects
*for stocks and the house, incorporating pension wealth
*It's used on the paper on wealth counterfactuals 
*for the verification for the Review of Economics and Statistics
*The program is run on Stata 11, using the version 9.2 option
*June 2011 - Dimitris Christelis

version 9.2
cap clear all
macro drop _all
clear
set mem 500m
set more off

*Paths
do "c:/foldp/foldp.do"
global projf "${foldp}/Wealth Counterfactuals/Verification"

*Tolerance
sca tolr=1e-6

*Weights
global wgt wgth

*ATTENTION: The base country (Germany in our case) has to be put first in the sequence
global eureg1 de se dk nl be fr ch at it es gr
global eureg2 Germany Sweden Denmark Netherlands Belgium France Switzerland Austria ///
   Italy Spain Greece
*Country numbers: they have to follow exactly the same sequence as the country names
global eureg3 3 1 2 4 5 6 7 8 9 10 11   
local nas: word count $eureg1

global assetl1 rfin home /*ownb secdebt liab tdebt*/
global assetl2 hrfin hhome /*hownb hsecdebt hliab htdebt*/
local nar: word count $assetl1

*The loop for the asset must come before the loop for the countries
 
forvalues ll=1/`nar' {

   *Globals of asset
   global as: word `ll' of $assetl1
   global hhas: word `ll' of $assetl2

   forvalues kk=1/`nas' {

   *Globals of country
   global re: word `kk' of $eureg1
   global fnre: word `kk' of $eureg2
   global cn: word `kk' of $eureg3

   *Logging
   cap log close
   log using "${projf}/Results/SHARE/${re}/penwealth/${as}/penwealth_${re}_${as}.log", replace


   *Loop for inclusion of pension wealth   
   forvalues www=1/3 {   
    
      if `www'==1 { 
         global pw nopw
         global aw 0
      }
      if `www'==2 { 
         global pw pw0
         global aw hpw0
      }
      if `www'==3 { 
         global pw pw1
         global aw hpw1
      }

   global pdate "$S_DATE" 
   global ptime "$S_TIME"
   display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"

   *Opening HRS data to generate quartiles
   use "${projf}/Data/HRS/HRS_DG.dta", clear

   *Checking the sample weights
   assert abs(wgth-wgtach)<1e-6
   drop wgtach

   *Keeping only the pension wealth sample
   qui keep if pwsample==1

   *Keeping only heads
   keep if head==1 & wgth>0 & wgth!=.

   *Putting pension wealth with a 1% real growth rate of pension benefits
   *equal to the one with 0% growth, so that it can be compared with the 
   *European countries
   qui gen double hpw1 = hpw0


   *Total debt
   qui gen double htdebtv = hmortv+hhmlov+hliabv
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv + hhmlov 
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v+${aw}
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v+${aw}
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Creating wealth quartiles 
   qui su nw if head==1 [aw=wgth], detail
   sca p25_w_us = r(p25)
   sca p50_w_us = r(p50)
   sca p75_w_us = r(p75)

   *Creating income quartiles
   qui su inc if head==1 [aw=wgth], detail
   sca p25_i_us = r(p25)
   sca p50_i_us = r(p50)
   sca p75_i_us = r(p75)


   *Reading the parameters, variable lists etc.
   do "${projf}/Programs/Parameters.do"

   *global nsim=2
   *global finsim=5*${nsim}

   *Reading the SHARE data
   use "${projf}/Data/SHARE/SHARE_DG.dta", clear

   *Keeping only the country of interest
   keep if country==${cn} & ${wgt}>0 & ${wgt}!=.

   *Keeping only the pension wealth sample
   qui keep if pwsample==1

   *Keeping only heads
   keep if head==1

   *Total debt
   qui gen double htdebtv = hmortv+hliabv /* +hhmlov*/
   qui gen byte htdebto = (htdebtv>tolr)

   *Secured debt, except on 2nd house
   qui ge double hsecdebtv = hmortv /*+ hhmlov + hmort2v*/
   qui ge double hsecdebto = (hsecdebtv>tolr)

   *Ownership dummy
   *qui gen byte ${hhas}o=(${hhas}v>tolr)
   
   *Rescaling age
   qui replace age=age/sc_age
   qui gen double agesq=age^2

   *Creating net worth and income quartiles
   *Net worth, equal total minus the asset
   if "${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb" {
      qui gen double nw=hnetwv-${hhas}v+${aw}
   }
   *For debts, we add back instead of substracting
   if ! ("${as}"=="rfin" | "${as}"=="home" | "${as}"=="ownb") {
      qui gen double nw=hnetwv+${hhas}v+${aw}
   }
   
   *Income, equal to total minus capital income
   qui gen double inc=hgtincv-hcpincv

   *Bad health dummy
   qui gen byte badhs=(srhealtha>=4)

   *College dummy
   qui gen byte college=(edu==3)
   
   *High school dummy
   qui gen byte highs=(edu==2)

   *Couple dummy
   qui gen byte couple=(mstatus==1)
   
   *Widow dummy
   qui gen byte widow=(mstatus==3)
   
   *Never married
   qui gen byte nevmar=(mstatus==4)
      
   *Creating wealth quartiles for financial and real assets, both with respect
   *to US and DE
   if ${cn}==3 {
      qui su nw if head==1 [aw=wgth], detail
      sca p25_w_de = r(p25)
      sca p50_w_de = r(p50)
      sca p75_w_de = r(p75)
   }
   qui gen double wqnw1_qus=(nw<=scalar(p25_w_us))
   qui gen double wqnw2_qus=(nw>scalar(p25_w_us) & nw<=scalar(p50_w_us))
   qui gen double wqnw3_qus=(nw>scalar(p50_w_us) & nw<=scalar(p75_w_us))
   qui gen double wqnw4_qus=(nw>scalar(p75_w_us))
   qui gen double wqnw1_qde=(nw<=scalar(p25_w_de))
   qui gen double wqnw2_qde=(nw>scalar(p25_w_de) & nw<=scalar(p50_w_de))
   qui gen double wqnw3_qde=(nw>scalar(p50_w_de) & nw<=scalar(p75_w_de))
   qui gen double wqnw4_qde=(nw>scalar(p75_w_de))
   
   if ${cn}==3 {   
      qui su inc if head==1 [aw=wgth], detail
      sca p25_i_de = r(p25)
      sca p50_i_de = r(p50)
      sca p75_i_de = r(p75)
   }
   qui gen double wqinc1_qus=(inc<=scalar(p25_i_us))
   qui gen double wqinc2_qus=(inc>scalar(p25_i_us) & inc<=scalar(p50_i_us))
   qui gen double wqinc3_qus=(inc>scalar(p50_i_us) & inc<=scalar(p75_i_us))
   qui gen double wqinc4_qus=(inc>scalar(p75_i_us))
   qui gen double wqinc1_qde=(inc<=scalar(p25_i_de))
   qui gen double wqinc2_qde=(inc>scalar(p25_i_de) & inc<=scalar(p50_i_de))
   qui gen double wqinc3_qde=(inc>scalar(p50_i_de) & inc<=scalar(p75_i_de))
   qui gen double wqinc4_qde=(inc>scalar(p75_i_de))


   *Saving the dataset
   qui save "${projf}/Results/SHARE/${re}/penwealth/${as}/temp.dta", replace

   *Running 2 regressions, using both definitions of quartiles
   foreach cb in qus /*qde*/ {


      *Bootstrap program

      *Creating a matrix with number of columsn equal to 1 + # bootstrap runs, storing
      *the main sample estimates in the first line, together with the proportion of
      *hhds owning the asset
      qui su ${hhas}o if head==1 [aw=$wgt], meanonly
      sca ${as}p=r(mean)
      mat ${as}coeffall=${as}p
      local nch=$nc+1
      mat colnames ${as}coeffall = _b`nch'   
      qui save "${projf}/Results/SHARE/${re}/penwealth/${as}/temp2.dta", replace
  
      *Subsequent runs, beginning of the bootstrap loop
      forvalues i=1/$finsim {

         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
           
         use "${projf}/Results/SHARE/${re}/penwealth/${as}/temp2.dta", clear
         *Fixing the seed before each resampling
         local sd=1234+30*`i'+`kk'
         set seed `sd'
         bsample
               
            qui su ${hhas}o if head==1 [aw=$wgt], meanonly
            sca ${as}p=r(mean)
            mat ${as}coeff`sr'=${as}p
            mat ${as}coeffall=(${as}coeffall \ ${as}coeff`sr')

            forvalues lll=1/10 {
               if `sr'==$nsim*(`lll'/10) {
                  noisily di in yellow `sr' in yellow "*" _continue 
               }
            }

            mat drop ${as}coeff`sr'

         
               
         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         if `sr'==$nsim {
            di in yellow "End of ${nsim} valid boostrap runs, total runs equal to `i'"
            continue, break
         }
         
      }  /* End of loop over bootstrap runs */

      *Saving the matrix of sample proportions
      drop _all
      svmat double ${as}coeffall
      qui save "${projf}/Results/SHARE/${re}/penwealth/${as}/boot_predpr_${re}_${as}_${pw}_`cb'.dta", replace

        
      ***********************************************************
      *Calculation of marginal effects and associated magnitudes*
      ***********************************************************

      global nc: word count ${basiccovs_`cb'} _cons
      
      *Reading the dataset
      drop _all
      use "${projf}/Results/SHARE/${re}/penwealth/${as}/temp.dta", clear
      sca nobs_`cb'=_N

      *Keeping only the pension wealth sample
      qui keep if pwsample==1

      *Incrementing the continuous variables by the appropriate magnitudes
      *to compute marginal effects and elasticities
      *Variables that are not used in elasticity form
      foreach y in $conlist_meff {
         qui gen double `y'_pl=`y'+(ch_`y'/sc_`y')
      }   

      *Merging with the sample proportions
      merge using "${projf}/Results/SHARE/${re}/penwealth/${as}/boot_predpr_${re}_${as}_${pw}_`cb'.dta"
      drop _merge
      
      *Renaming the coefficient of the constant 
      *set trace on
      rename ${as}coeffall pr_${re}_${as}_${pw}_`cb'_own

   
      *Merging with the draws from the US
      if "`cb'"=="qus" {
         merge using "${projf}/Results/HRS/Country/penwealth/${as}/boot_drawcoeff_usall_${as}_${pw}.dta"   
      	 drop _merge   

      	 *Putting the values of the drawn coefficients equal to the dummy value from missing
      	 *in rows lower than the number of bootstraps, so that there are no missing values
      	 local nch = $nc + 1
      	 forvalues h=1/`nch' {   
      	    qui replace _b`h'=$dumval if _n>${nsim}+1  
      	 }      
   
      	 *Renaming the coefficients  
      	 local ncl = $nc - 1
      	 forvalues h=1/`ncl' {   
      	    local g: word `h' of ${basiccovs_qus}
      	    rename _b`h' b_usall_`g'   
      	 }   

      	 *Renaming the coefficient of the constant 
      	 *set trace on
      	 rename _b${nc} b_usall_cons
   
      	 local nch = $nc + 1
      	 rename _b`nch' pr_usall_${as}_${pw}_own
      }
      
      *Merging with the coefficients and the probabilities from Germany
      if "`cb'"=="qde" {
         if "${re}"~="de" {

      	    *Coefficients
      	    merge using "${projf}/Results/SHARE/de/${as}/boot_drawcoeff_de_${as}_${pw}_qde.dta"   
      	    drop _merge   

      	    *Putting the values of the drawn coefficients equal to the dummy value from missing
      	    *in rows lower than the number of bootstraps, so that there are no missing values
      	    local nch = $nc + 1
      	    forvalues h=1/`nch' {   
      	       qui replace _b`h'=$dumval if _n>${nsim}+1  
      	    }      
   
      	    *Renaming the coefficients  
      	    local ncl = $nc - 1
      	    forvalues h=1/`ncl' {   
      	       local g: word `h' of ${basiccovs_qde}
      	       rename _b`h' b_de_`g'_qde   
      	    }   

      	    *Renaming the coefficient of the constant 
      	    *set trace on
      	    rename _b${nc} b_de_cons_qde
   
      	    local nch = $nc + 1
      	    rename _b`nch' pr_de_${as}_${pw}_qde_own
      	 }         
      }
      
      *Loop over bootstrap coefficients
      forvalues i=1/$finsim { 

         *Defining an index of valid bootstrap runs that will determine whether the 
         *loop finishes or not
         if `i'==1 {
            local sr = 1
         }
         *Augmenting the index that tracks valid bootstrap runs by one
         if `i'>1 {
            local sr = `sr' + 1
         }
         *If there were any condition on the coefficients that was violated by the current
         *draw then the index of valid bootstrap runs would be decreased by one at this point
         *in the program. For the pooled probit there is no such condition and thus the
         *index of valid bootstrap runs (sr) is equal to the on of total bootstrap runs (i)
   
      
         *Creating a counterfactual prediction using the coefficients from the US
         if "`cb'"=="qus" {
            qui gen double xb0_usall = b_usall_cons[`i']
            foreach y in $basiccovs_qus {
                qui replace xb0_usall = xb0_usall + b_usall_`y'[`i']*`y' 
            }
	 }

         *Creating a counterfactual prediction using the coefficients from Germany
         if "`cb'"=="qde" {
            if "${re}"~="de" {         
               qui gen double xb0_de_qde = b_de_cons_qde[`i']
               foreach y in $basiccovs_qde {
                  qui replace xb0_de_qde = xb0_de_qde + b_de_`y'_qde[`i']*`y' 
               }
            }            
         }
      
         *******************************************************************
         *Calculating the average of probabilities across individual obs
         if "`cb'"=="qus" {   
            qui gen double dp_usall=normal(xb0_usall)
            qui su dp_usall [aw=$wgt], meanonly
            sca r0m=r(mean)
         }
         
         *Probability using Germany coefficients
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               qui gen double dp_de_qde=normal(xb0_de_qde)
               qui su dp_de_qde [aw=$wgt], meanonly
               sca r2m_qde=r(mean)
            }
         }

         if `i'==1 {
             
            if "`cb'"=="qus" {
               *Counterfactual probability using US coefficients
               qui gen double pr_${re}_${as}_${pw}_qus_cus=r0m
             
               *Difference between US ownership proportion and
               *own ownership proportion
               qui gen double dpus_${re}_${as}_${pw}_qus_tot=pr_usall_${as}_${pw}_own - ///
                  pr_${re}_${as}_${pw}_qus_own
             
               *Difference between US and own probability calculated using
               *US coefficients (characteristics effect)
               qui gen double dpus_${re}_${as}_${pw}_qus_char=pr_usall_${as}_${pw}_own[1] - ///   
                  pr_${re}_${as}_${pw}_qus_cus[1]
             
               *Difference between own probability calculated using
               *US coefficients and own probability (coefficient effect)
               qui gen double dpus_${re}_${as}_${pw}_qus_coeff=pr_${re}_${as}_${pw}_qus_cus[1] - ///   
                  pr_${re}_${as}_${pw}_qus_own[1]
            } 

            if "`cb'"=="qde" {
               if "${re}"~="de" {
            
             	  *Counterfactual probability using German coefficients
             	  qui gen double pr_${re}_${as}_${pw}_qde_cde=r2m_qde
             
             	  *Difference between mean German predicted probability and
             	  *mean own predicted probability
             	  qui gen double dpde_${re}_${as}_${pw}_qde_tot=pr_de_${as}_${pw}_qde_own - ///
             	     pr_${re}_${as}_${pw}_qde_own
             
             	  *Difference between German and own probability calculated using
             	  *German coefficients (characteristics effect)
             	  qui gen double dpde_${re}_${as}_${pw}_qde_char = pr_de_${as}_${pw}_qde_own[1] - ///
             	     pr_${re}_${as}_${pw}_qde_cde[1]  
             
             	  *Difference between own probability calculated using
             	  *German coefficients and own probability (coefficient effect)
             	  qui gen double dpde_${re}_${as}_${pw}_qde_coeff = pr_${re}_${as}_${pw}_qde_cde[1] - ///
             	    pr_${re}_${as}_${pw}_qde_own[1]
               }
            }
         }
         
         *Putting the probability at the line indexed by
         *the valid bootstrap run (sr), and not the general
         *run (i)
         
         if `i'>1 {
            if "`cb'"=="qus" {
               
               *Counterfactual probability using US coefficients   
               qui replace pr_${re}_${as}_${pw}_qus_cus=r0m in `sr'   
               assert pr_${re}_${as}_${pw}_qus_cus<$dumval if _n==`sr'   

               *Difference between US and own probability calculated using   
               *US coefficients (characteristics effect)   
               qui replace dpus_${re}_${as}_${pw}_qus_char=pr_usall_${as}_${pw}_own[`sr'] - ///   
                  pr_${re}_${as}_${pw}_qus_cus[`sr'] in `sr'   

               *Difference between own probability calculated using   
               *US coefficients and own probability (coefficient effect)   
               qui replace dpus_${re}_${as}_${pw}_qus_coeff=pr_${re}_${as}_${pw}_qus_cus[`sr'] - ///   
                  pr_${re}_${as}_${pw}_qus_own[`sr'] in `sr'   
            }

            if "`cb'"=="qde" {
               if "${re}"~="de" {
                  *Counterfactual probability using German coefficients
               	  qui replace pr_${re}_${as}_${pw}_qde_cde=r2m_qde in `sr'
               	  assert pr_${re}_${as}_${pw}_qde_cde<$dumval if _n==`sr'
         
               	  *Difference between German and own probability calculated using
               	  *German coefficients (characteristics effect)
               	  qui replace dpde_${re}_${as}_${pw}_qde_char=pr_de_${as}_${pw}_qde_own[`sr'] - ///
               	     pr_${re}_${as}_${pw}_qde_cde[`sr'] in `sr'

               	  *Difference between own probability calculated using
               	  *German coefficients and own probability (coefficient effect)
               	  qui replace dpde_${re}_${as}_${pw}_qde_coeff=pr_${re}_${as}_${pw}_qde_cde[`sr'] - ///
               	     pr_${re}_${as}_${pw}_qde_own[`sr'] in `sr'
               } 
            }
         }
         
         if "`cb'"=="qus" {
            drop dp_usall xb0_usall
         }
         
         if "`cb'"=="qde" {
            if "${re}"~="de" {
               drop dp_de_qde xb0_de_qde
            }   
         }   
      
         forvalues lll=1/10 {
            if `sr'==$nsim*(`lll'/10) {
               di in yellow "Iteration `sr', ${as}"
            }
         }

         *If the number of valid bootstrap runs is equal to the predetermined
         *one then the loop over bootstrap runs stops
         local ncpp=$nsim+1
         if `sr'==`ncpp' {
            di in yellow "End of ${nsim} valid bootstrap runs, total runs equal to `i'"
            continue, break
         }
      
      }  /* End of loop over bootstrap runs */

      *Dropping obs beyond the line of the last valid bootstrap run
      drop if _n>$nsim+1
   
      *Global of marginal effects, putting first the probabilities of interest
      if "`cb'"=="qus" {
         local me pr_${re}_${as}_${pw}_qus_own pr_${re}_${as}_${pw}_qus_cus dpus_${re}_${as}_${pw}_qus_tot ///
            dpus_${re}_${as}_${pw}_qus_coeff dpus_${re}_${as}_${pw}_qus_char
      }
      if "`cb'"=="qde" {
         if "${re}"=="de" {
            local me pr_${re}_${as}_${pw}_qde_own
         } 
         if "${re}"~="de" {
            local me pr_${re}_${as}_${pw}_qde_own pr_${re}_${as}_${pw}_qde_cde dpde_${re}_${as}_${pw}_qde_tot ///
               dpde_${re}_${as}_${pw}_qde_coeff dpde_${re}_${as}_${pw}_qde_char
         } 
      }
      

      *Keeping only the marginal effects and percentage variation
      keep country `me' `perc' 

      *Calculating the variance of the generated average probabilities and marginal effects
      *One has to start from the 2nd line, since the first line gives the magnitudes
      *evaluated at the estimated coefficients
      local ncpp=${nsim}+1
      foreach vr in `me' `perc' {
         qui su `vr' in 2/`ncpp' 
         qui gen double se_`vr'=sqrt(r(Var))
      }

      *Saving the dataset with all the probabilities and marginal effects
      qui save "${projf}/Results/SHARE/${re}/penwealth/${as}/boot_alldraweff_${re}_${as}_${pw}_`cb'.dta", replace

      *Saving the dataset with only the mean magnitudes and their standard errors
      qui keep in 1
      qui save "${projf}/Results/SHARE/${re}/penwealth/${as}/boot_meaneff_${re}_${as}_${pw}_`cb'.dta", replace
 
      *Saving matrices of marginal effects and percentage variations as a text file

      foreach rs in $prfx  {
         *Number of elements
         local n_`rs': word count ``rs''
         matrix olout=J(`n_`rs'',6,.)

         forvalues k=1/`n_`rs'' {
            local h: word `k' of ``rs''
            matrix olout[`k',1] = `h'[1]
            matrix olout[`k',2] = se_`h'[1]
            matrix olout[`k',3] = olout[`k',1]/olout[`k',2]
            matrix olout[`k',4] = 2*ttail(nobs_`cb' - ${nc},abs(olout[`k',3]))
            matrix olout[`k',5] = olout[`k',1] - invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
            matrix olout[`k',6] = olout[`k',1] + invttail(nobs_`cb' - ${nc},sigl/2)*olout[`k',2]
         } 
    
         *Labels of matrix of magnitudes
         *Column labels
         local colmat `rs' se tstat pvalue ci_low ci_high
         matrix colnames olout = `colmat' 
         *Row labels
         matrix rownames olout = ``rs''
         *Displaying and saving the matrix of results
         mat list olout,format(%15.5f)
         mat2txt, mat(olout) ///
            saving("${projf}/Results/SHARE/${re}/penwealth/${as}/boot_mefftab_${re}_${as}_${pw}_`cb'.txt") replace format(%15.5f) ///
            title("Meff for ${as} in ${fnre}, using `cb' quartiles")

      }  /* End of loop over magnitudes */        

      erase "${projf}/Results/SHARE/${re}/penwealth/${as}/temp2.dta" 
   } /* End of loop over quartile specifications */
  
   erase "${projf}/Results/SHARE/${re}/penwealth/${as}/temp.dta"         

   }  /* End of loop over the kinds of pension wealth */

   }  /*End of loop over countries */    
}  /*End of loop over assets */    


display "log file for ${as} in ${fnre} opened on ${pdate} at ${ptime}"
display "log file for ${as} in ${fnre} closed on $S_DATE at $S_TIME"
cap log close
   
   
  


